# -*- coding: utf-8 -*-
"""
Created on Mon Sep 16 15:50:24 2013

@author: jkwong
"""

#PBAR_ZspecBasic_ExamineAirSection.py

import PBAR_Zspec
import numpy as np
import os
import PBAR_Cargo


# Load the data

acquisitionTime = np.array([1/60., 1/60.])
airSection = np.array([0, 175])

filenameList = []
filenameList.append('4939-FDFC-All2157.npy')
filenameList.append('4955-FDFC-All2154.npy')

fullFilenameList = []
fullFilenameList.append(r'C:\Users\jkwong\Documents\Work\PBAR\data3\4939\4939-FDFC-All2157.npy')
fullFilenameList.append(r'C:\Users\jkwong\Documents\Work\PBAR\data3\4955\4955-FDFC-All2154.npy')

energyList = []
dat = []
for (index, fullFilename) in enumerate(fullFilenameList):
    (a, b) = PBAR_Zspec.ReadZspecBasicScanNumpy(fullFilename)
    energyList.append(a)
    dat.append(b)
del a
del b

# load the corresponding zspec closed collimator calibration data
acquisitionTimeCC = np.array([300., 300.])
filenameCCList = []
filenameCCList.append('eb16.csv')
filenameCCList.append('eb17.csv')
baseDir = r'C:\Users\jkwong\Documents\Work\PBAR\data';

fullFilenameCCList = np.array([os.path.join(baseDir, f) for f in filenameCCList])

datCC = []
datCC = PBAR_Zspec.ReadZspec(fullFilenameCCList)

# Calculate the discriminants (will move this to PBAR_Zspec package)

discrim = []

for (index, d) in enumerate(dat):
    temp = {}
    energy = energyList[index]
    
    temp['count'] = d.sum(2)
    
    energyMatrix = np.tile(energy, (d.shape[0], d.shape[1], 1))
    temp['binMean'] = (d * energyMatrix).sum(2) / d.sum(2)
    
    binMeanMatrix = np.tile(temp['binMean'], (d.shape[2], 1, 1))
    temp['binSTD'] = np.sqrt( (( (energyMatrix - binMeanMatrix)**2) * d ).astype(float) / d.sum(2).astype(float)  )

    temp['binSTD_binMean'] = temp['binSTD'] / temp['binMean']

# show the images

discrimName = 'binMean'



#plt.imshow(intensity, interpolation = 'nearest', aspect='auto', cmap = cm.Greys_r)

# subplots of spectra of all good detectors


# Plot the waveforms at multiple times bins, SUBPLOTS

# index = 1
#index = 6
peakIndex = 0

plotFits = 1
filterSpectra = 0
filterWidth = 1
normalize = 0
correctGain = 0
plotRate = 1

subplotx = 4
subploty = 2

# cycle through detectors
figNum = 0
for detNum in np.arange(100,120):
    # Generate a new figure window if filled.
    if ((figNum % (subplotx * subploty)) == 0):
        f, ax = plt.subplots(subploty, subplotx, sharex='col', sharey='row')
    # Calculate subplot window index
    ii = (figNum % (subplotx*subploty)) /subplotx
    jj = (figNum % (subplotx*subploty)) % subplotx

    
    for i in xrange(0, 2):
        energy = energyList[i]
        if correctGain:
            x = energy * gainShift[i,detNum]
        else:
            x = energy
            
        if plotRate:
            ax[ii][jj].plot(x, \
                dat[i][airSection[0]:airSection[1], detNum,:].sum(0)/acquisitionTime[i]/(airSection[1] - airSection[0]) + 1e-10, \
                label = '%s' %filenameList[i])
        else:
            ax[ii][jj].plot(x, \
                dat[i][airSection[0]:airSection[1]][detNum,:].sum(0)/acquisitionTime[i]/(airSection[1] - airSection[0]) + 1e-10, \
                label = '%s' %filenameList[i])

    for i in xrange(0, 2):
        if correctGain:
            x = np.arange(256).astype(float) * gainShift[i,detNum]
        else:
            x = np.arange(256).astype(float)
        if plotRate:
            ax[ii][jj].plot(x, \
                datCC[i][:,detNum]/acquisitionTimeCC[i] + 1e-10, \
                label = '%s' %filenameCCList[i])
        else:
            ax[ii][jj].plot(x, \
                datCC[i][:,detNum]/acquisitionTimeCC[i] + 1e-10, \
                label = '%s' %filenameCCList[i])

#        ax[ii][jj].plot(x, dat[i][:,detNum]/acquisitionTime[i]/pulseRate[i] + 1e-10, label = '%s' %filenameList[i])
    
    ax[ii][jj].grid()
    if (ii == 0) & (jj == 0):
        ax[ii][jj].legend()
    #ax[ii][jj].set_xlim((100, 200))
#    ax[ii][jj].set_xlim((60, 160))
    
    ax[ii][jj].set_yscale('log')
    if plotRate:
        ax[ii][jj].set_ylabel('Rate')
#        ax[ii][jj].set_ylim((1e-2, 1e1))
    else:
        ax[ii][jj].set_ylabel('Count')
#        ax[ii][jj].set_ylim((1e-2, 1e6))
    ax[ii][jj].set_xlabel('Bin')
    ax[ii][jj].set_title('Det# %d' %(detNum+1))
    figNum += 1


# plot of binMean
